1. Introduction to Finite Element Method 


Instructional Objectives of Chapter 1 


1. Provide the background of the finite element method (FEM). 


2. Explain how the finite element method is used in practice. 


There are quite a few engineering problems for which an exact solution does not exist, 
and to find its closed-form may be cumbersome to implement or too expensive. Some 
reasons that could make the exact solution inadequate include the complicated nature of 
differential equations, geometry, and the challenges arising from the boundary and/or initial 
conditions. Numerical approximations help overcome these challenges. The following are 
some of the spatial/temporal discretization methods used in computational mechanics, and 
they all consist of discretizing the system using governing partial differential equations 
representing the physical system. 


1. Finite Difference Method (FDM): The finite difference method is a direct method 
for approximating partial differential equations on a discrete grid by approximating 
derivatives of the unknown quantities grid by linear differences. This results in a 
set of simultaneous linear equations. Although these methods are easy to grasp and 
implement, they are burdensome when applied to problems involving complicated 
geometries or boundary conditions. This method remains the top preference for 
many fluid mechanics problems. 
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2. Finite Element Method (FEM): The finite element method uses integral formulations 
rather than difference equations to create a system of algebraic equations to approxi- 
mate solutions to the differential equation. In this method, the domain is subdivided 
into small subdomains, and each subdomain uses polynomials to approximate the 
solution. 


3. Boundary Element Method (BEM): The boundary element method is a numerical 
computational method where the solution of linear partial differential equations 
involves the boundary integral formulation. The boundary element method attempts 
to use the given boundary conditions to fit boundary values into the integral equation 
rather than values throughout the space defined by a partial differential equation. 


4. Finite Volume Method (FVM): The finite volume method is a numerical method 
for solving partial differential equations where the calculations average the values 
of the conserved variables across the volume. One advantage of the finite volume 
method over finite difference methods is that it does not require a structured mesh. 
The finite volume method is preferable to other methods because the application of 
the boundary conditions can be noninvasive. 


5. Meshfree Method (MM): Meshfree methods is a method where the development (J. Awrejcewicz, 2011) 
of the system algebraic equations for the whole domain of the problem does not 
use a predefined mesh for the domain discretization. This technique uses a set of 
scattered nodes, called field nodes, to establish the problem domain and boundaries, 
which do not require any prior information on the relationship between the nodes 
for the interpolation or approximation of the unknown functions of field variables. 


This book focuses on the finite element method, and this chapter introduces the 
finite element method as the numerical method to approximate solutions to engineering 
problems. 


1.1 Finite Element Method 


The design of many modern engineering products heavily uses the Finite Element Method 
(FEM). These methods can be used along with physical testing to determine the product's 
capability to withstand the environment intended. Testing is extremely costly, and in some 
cases, it may not be possible to achieve target test loads in all critical areas of the structure. 
FEM can predict the system's behavior for design validation, anomaly investigations, 
improve test coverage in test programs, and analyze regions of the structure where it may 
not be possible to achieve the target loads. 


FEM has its origins in mathematical principles. Although Courant provided the 
earlier mathematical treatment of the method in 1943, the application of the method 
did not take place until 1968. This first application was for electromagnetic problems, 
ever since expanded to diverse fields such as waveguide problems, electric machines, 
semiconductor devices, microstrips, biological bodies, and structural applications. 
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The finite difference method (FDM) is conceptually more straightforward and easier 
to implement into a computer program when compared to the FEM. However, FEM is a 
powerful and versatile numerical technique that can handle complicated geometries, inho- 
mogeneous media, nonlinearities, transient problems, and discontinuities. The systematic 
nature of the method makes it possible to construct general-purpose computer programs 
for a wide range of problems. The computer programs can easily be adapted to different 
fields with little or no modification. 


What is FEM? 


FEM is a numerical technique that approximates the solution to partial differential equa- 
tions. Finite element analysis (FEA) is the application of FEM to approximate solutions to 
real problems involving complicated geometry, complicated materials properties, and/or 
complicated boundary or initial conditions. FEM is used in most engineering design and 
analysis from the early stages of the design process all the way into the production phase. 
Table 1.1 shows some typical engineering applications that use FEM. 


In-depth knowledge of the finite element theory is necessary to develop finite element 
codes, and basic knowledge is required to use the FEA software properly. The knowledge 
gained dramatically enhances the engineers’ capacity to find modeling errors and to 
interpret results correctly. Skills in FEA can enhance an engineer’s career to aid them in 
the design process, component-level analysis, product development, or applied research. 


FEA Background 


Some argue that variational calculus formed the basis for the FEM. In 1696, Bernoulli 
developed the foundations of variational calculus. However, in 1744, Euler gave a general 
solution to variational problems for a differential equation. Eleven years later, Lagrange 
proposed Euler an elegant justification for this equation. The prodigious contributions of 
Euler concerning the analytic and numerical solutions to differential equations form the 
foundations of FEM. 


The origins of FEM philosophy dates back to the late 1800s and early 1900s when 


Table 1.1: Typical FEM applications 


Application Primary variable 

Heat transfer T (x,y,z) 
Solid mechanics | u(x,y),v(x, y) 
Fluid mechanics v(x,y,z) 

Electric fields E (r) 
Magnetic fields B(r) 
Acoustic fields p,u 
Coupled physics — 
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some researchers approximated and modeled elastic continua using discrete equivalent 
elastic bars. Rayleigh developed a method to predict the first natural frequency of simple 
structures. Later, Ritz proposed numerical analysis using variational calculus, expanding 
the method into the well-known Rayleigh-Ritz method to predict structural stress and 
displacement behaviors. However, the choice for the assumed shape, critical to the 
accuracy of the results, proved to be too difficult for complicated geometries. In 1963, 
mathematicians recognized the Rayleigh-Ritz technique, leading to the development of 
more robust FEA algorithms. 


As used today, FEM represents the convergence of three areas: matrix structural  (Felippa, 2007) 
analysis (MSA), variational approximation theory, and the digital computer. These did 
not appear simultaneously but came together in the early 1950s. MSA came on the — (Duncan and Collar, 1934) 
scene in the mid-1930s when desk calculators became popular. Before the age of digital (Duncan and Collar, 1935) 
computers, mathematicians proposed variational approximation schemes similar to those 
used in FEM. 


In 1941, Hrenikoff developed the “framework method” by formulating two-dimensional 

structural problems as an assembly of bars and beams. Later in 1943, R. Courant applied 

the Ritz method to obtain “piecewise approximations” to problems involving equilibrium 

and vibrations. One of Courant’s most significant contributions is the approximation (Courant, 1943) 
of the torsional rigidity of a hollow shaft by dividing the cross-section into triangles and 

linearly interpolating the stress function over each triangle. Further development of these 

ideas continued through the 1940s and 1950s. Its popularity grew in the 1950s when 

engineers began to use computers to solve structural problems. 


FEM was further developed in the 1950s by Boeing and Bell Aerospace in the United 
States and Rolls Royce in the United Kingdom. The first publications in the field were 
by Turner and his team. M. J. Turner, a Boeing employee at the time, generalized and (Turner et al., 1956) 
perfected the Direct Stiffness Method and got Boeing to commit resources to it, while 
other aerospace companies kept using the Force Method. During 1952-53, he oversaw 
the development of the first continuum-based finite elements. In 1956, Turner derived 
stiffness matrices for structures such as trusses, beams, among others. 


In addition to Turner, there were other major contributors to the current practice. B. 
M. Irons invented isoparametric models, shape functions, the patch test, and frontal solvers. 
R. Clough, a professor at the University of California at Berkely, led further maturity 
of the FEM. R. J. Melosh recognized that the Rayleigh-Ritz could systematically assist 
in the derivation of stiffness elements using variational principles. Along with Argyris’ 
program, they prototyped the first-generation FEM from 1950 through 1962. E. L. Wilson, — (Argyris and Kelsey, 1960) 
professor at Berkeley, developed one of the first FEA programs, and it became popular in 
the 1960s as it was a freeware program. 


These pioneers spent part of the career in the aerospace industry, and they were 
structural engineers with expertise in classical mechanics and applied mathematics. Since 
the advancement in digital computation is a crucial foundation for FEA, the developers 
of FEM interacted closely with the aircraft industry, as they had the resources to tackle 
these computationally challenging problems. Since large aerospace companies used FEM, 
the focus at the time was on thin structures built up with bars, ribs, spars, stiffeners, and 
panels. (Felippa, 2007) 
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In 1960, following the lead of Clough, researchers and engineers formally started to 
use the term finite element. In 1963, finite element analysis became more popular, leading 
to its expansion and in-depth research. This led to its application beyond structures; to 
fields such as heat transfer and fluid mechanics problems. In 1964, Clough convinced 
Olek Zienkiewicz, originally an expert in finite difference methods, to attempt the FEM to 
solve problems. Zienkiewicz and Cheung (1967) wrote the first book on the subject and 
organized another important Civil Engineering research group at the University of Wales 
at Swansea. 


By the early 1970s, finite element analysis had established itself as a general numer- 
ical technique for solving any system of partial differential equations. With the advent 
of powerful computers, the method kept expanding its horizons to many industries such 
as aerospace, automotive, defense, and nuclear industries. In the early 2020s, typical 
simulations involved as many as tens of millions of degrees of freedom. 


1.1.3 Morphing from MSA to FEM 


Figure 1.1 shows that Matrix Structural Analysis (MSA) is the ancestor of FEM. MSA 

morphed from the pre-computer era into the first programmable computers through wobbly 

gyrations. MSA and FEM stand on three legs: mathematical models, matrix formulation (Felippa, 2001) 
of the discrete equations, and computing tools to do the numerical work. Of the three legs, 

the latter is the one that has undergone the most dramatic changes. The “human computers" 

of the 1930s and '40s morphed into programmable analog and digital computers. 


The mid-1960s settled the FEM configuration shown on the right in Fig. 1.1. The 
choice of primary unknowns has traditionally classified matrix formulations for MSA 
and FEM. In the Displacement Method (DM), the primary unknowns are the physical or 
generalized displacements. In the Classical Force Method (CFM), these are amplitudes 


MATRIX 
STRUCTURAL 
ANALYSIS 


Matrix Human 
Formulations Computers 


Figure 1.1: Morphing of the pre-computer MSA (before 1950) into the present FEM. On 
the left “human computer” means computations under direct human control, possibly with 
analog devices (slide rule) or digital devices (desk calculator). (Felippa, 2007) 
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of redundant force (or stress) patterns. Additional methods involve the combinations of 
displacements, forces, and/or deformations as primary unknowns. Nowadays, the term 
Stiffness Method replaced DM, and Flexibility Method replaced CFM. Sometimes the 
Stiffness Method is also called Direct Stiffness Method (DSM). 


1.1.4 COTS FEA Software 


All commercial off-the-shelf or commercially available off-the-shelf (COTS) FEA software 
have their unique purposes and advantages. It is hard to choose between one software 
over the other. Typically, it is the application that dictates the choice. Among the existing 
COTS FEA Software, the following are the most common ones: 


1. In the 1960s, Westinghouse Electric Corp engineer John Swanson developed soft- 
ware for nuclear reactor analysis, known as ANSYS. ANSYS is a large-scale 
general-purpose finite element analysis software that integrates structure, heat, fluid, 
electric field, magnetic field, and sound field analysis. ANSYS can provide the rich- 
ness of functionality across a broad range of disciplines, whether explicit, structural, 
fluid, thermal, or electromagnetics. Allows for truly coupling multiple physics in 
a single simulation and adapts to different CAD, PLM, in-house codes, and other 
point solutions typically comprise the overall design and development process. 


2. In 1965, NASA funded a project lead by Dick MacNeal to develop NASTRAN. 
Today, MSC-NASTRAN is a comprehensive and widely used general structural FEA 
software. It can effectively solve problems involving strength, stiffness, buckling, 
modal, dynamics, thermodynamics, nonlinearity, (noise) acoustics, fluid-structure 
coupling, aeroelasticity, supercell, inertial release, and structural optimization of 
various large and complex structures. 


3. In 1978, HKS developed ABAQUS. ABAQUS became, and still is, very popular 
among the research community because it allowed its users to develop new material 
models and elements. The Abaqus Unified FEA (previously called ABAQUS) can 
include full vehicle loads, dynamic vibration, multibody systems, impact/crash, 
nonlinear static, thermal coupling, and acoustic-structural coupling while using a 
standard model data structure and integrated solver technology. 


4. In 1989, John Hallquist developed LS-DYNA for crashworthiness, sheet metal 
forming, and prototype simulations. LS-DYNA is an advanced general-purpose 
multiphysics simulation software package developed by Livermore Software Tech- 
nology Corporation (LSTC) and acquired by ANSYS in 2019. The software can 
calculate many complex, real-world problems. Its core competency lies in highly 
nonlinear transient dynamic FEA using explicit time integration. 
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(a) (b) (c) 


Figure 1.2: Figure shows the discretization of: (a) a car (Mit, 2021), (b) a human joint 
(Beidokhti et al., 2017), and (c) a F-16 aircraft (Logston, 2021). 


1.2 The FEM Analysis Process 


So far, a background of the FEM was provided. The intent here is to provide an overview 
of how the method works. FEM is a systematic process consisting of dividing a domain 
into subdomains, called elements. Each element is analyzed individually. Figure 1.2 
shows the discretization of three real structures. Generally, it is possible to divide the 
structure into smaller components or geometrical shapes. Mathematically, the discretized 
continuous structure is now a set of finite subdomains or elements, leading to the name 
finite element analysis. Note how each finite element interconnects at more than one node. 
Boundary lines or planes define the subdomain. Instead of analyzing the entire system 
at once, the focus shifts to analyzing each element using the governing equations of the 
problem. After completing the analysis of each element and formulating it only in terms 
of unknown quantities at the nodes, then it is possible to find its solution by solving the 
resulting system of equations by assembling the elements. This encapsulates the three 
overarching steps in the Finite Element Method: preprocessing (discretization strategy), 
processing (solution execution), and postprocessing (dissemination through data and visual 
representation of the results). 


1.2.1 Illustration of Discretization Concepts 
(Carlos, 2007) 


In order to understand how discretization concepts work, consider the circle of radius Example Courtesy of Dr. Carlos Felippa, Univer- 
r, as in Fig. 1.3(a). It is possible to use an ancient problem to illustrate the underlying sity of Colorado at Boulder. 

concept: find the perimeter p of a circle with diameter d. Since L = zd, the problem 

is equivalent to finding the numerical value for 7. Now, inscribe a regular polygon of n 

sides, where n = 8 in Fig. 1.3(b). Rename polygon sides as elements and vertices as nodes. 

Label the nodes with integers 1,...,8. Extract a typical element, say that joining nodes 

4-5, as in Fig. 1.3(c). This is an instance of the generic element i—j pictured in Fig. 1.3(d). 

The element length is Li; = 2rsin(z//n). Since all elements have the same length, the n* 
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(a) (b) (c) (d) 
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Figure 1.3: The “find 7” problem treated with FEM concepts: (a) continuum object, (b) 
a discrete approximation by inscribed regular polygons, (c) disconnected element, (d) 
generic element. (Felippa, 2007) 


polygon perimeter is Lj, = n Lij, whence the approximation to 7 is 


Ln s fT : 
Tn = d =nsin| —), (nisan even number) 


Table 1.2 summarizes the values of m, for n = 1,2,4,...,10 000 000, showing are 
fairly slow convergence to 7. The last value displays 15-place accuracy. As the number 
of terms keeps increasing, the solution approaches the exact solution, i.e., the discretized 
solution approaches a circle. This simple example shows some key ideas behind the 
mathematical discretization: 


1. Preprocessing: The circle is a source mathematical object that the polygons replace. 
These are discrete approximations to the circle. The sides, renamed as elements, are 
specified by their end nodes. 


Table 1.2: Rectification of Circle by Inscribed Polygons 


n Tn —nsin(z/n) Exact 7 to 16 places 
1 0.000000000000000 
2 2.000000000000000 
4 2.828427124746190 
8 3.061467458920718 
16 3.121445152258052 
32 3.136548490545939 
64 3.140331156954753 
128 3.141277250932773 

256 3.141513801144301 | 3.141592653589793 


10,000 3.141592601912665 | 3.141592653589793 
100,000 3.141592653073022 | 3.141592653589793 
10,000,000 | 3.141592653589741 | 3.141592653589793 
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2. Processing: Elements can be separated by disconnecting nodes, a process called 
disassembly in the FEM. Upon disassembly, a generic element can be defined, 
independently of the original circle, by the segment that connects two nodes i and 
j. The relevant element property: side length Lij, can be computed in the generic 
element independently of the others, a property called local support in the FEM. 
The target property: the polygon perimeter, is obtained by reconnecting n elements 
and adding up their length; the corresponding steps in the FEM being assembly and 
solution, respectively. 


3. Postprocessing: Table 1.2 shows the final results. 


There is nothing magical about the circle; the same technique is applicable to rectify any 

smooth plane curve. Martin and Carey (1973) used this example to show that finite (Martin and Carey, 1973) 
element ideas are traceable back to the Egyptian mathematicians from circa 1800 B.C. and 

Archimedes’ famous studies on circle rectification by 250 B.C. However, a comparison 

with the modern FEM may be a bit of a stretch because the example does not illustrate 

the concept of degrees of freedom, conjugate quantities, and local-global coordinates. It 

is also guilty of circular reasoning: the compact formula z = lim,_,.. n sin(z/n) uses 

the unknown z in the right-hand side. Despite these flaws, the example is useful in one 

respect: showing it is possible to replace one mathematical object for another. 


1.2.2 Systematic Approach to the Finite Element Method 


(Carlos, 2007) 


PHYSICAL 
SYSTEM 


VERIFICATION 
SIMULATION ERROR 


VALIDATION 
SIMULATION ERROR: MODELING & SOLUTION ERROR 


Figure 1.4: The Physical FEM. The physical system (left) is the source of the simulation 
process. The ideal mathematical model (should one go to the trouble of constructing it) is 
inessential. (Felippa, 2007) 


Modeling complicated physical systems such as an aircraft subjected to aerodynamic 
loads requires a systematic approach to solving the problem numerically. The formulation 
of the numerical techniques must be such that they can solve any complicated problem 
within a computer program. The idea is to use a model-based simulation to solve compli- 
cated physical systems, Figure 1.4. As with the circular example, the process consists of 
idealizing and discretizing the problem to produce results for the problem of interest. 
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Usually, the representation of the behavior of physical models is through mathematical 
models consisting of an ordinary or partial differential equation in space and time. The 
analytical solutions to these problems are impractical due to the complicated nature of 
the geometry and boundary conditions. Physical problems can be structural, heat transfer, 
fluids, and electrostatic charge-distribution problems. Simplifications to these physical 
problems have the purpose of reducing computational time. These idealizations are 
necessary to make the solution to these problems computationally tractable. 


1.2.3 Model Idealization 


Many engineering problems involve complicated geometries. In many cases, too dense 
of a finite element mesh may be required to capture the physical behavior of interest. 
Increasing the mesh density increases the number of degrees of freedom, making it 
computationally expensive. One alternative is to simplify the physical model using a 
combination of idealizations and boundary conditions to capture the real system’s behavior. 
Figure 1.5 shows two case scenarios of model idealization. The first case is the model of a 
bonded-joint of a wing spar. Here, the idealized physical model used in FEM results from 
modeling a portion of the web and making approximations to the boundary conditions 
results. The second case is the model of a centered-hole plate subject to a tensile axial 
load. Modeling only a fourth of the model with boundary conditions addressing symmetry 
results in the idealized physical model used in FEM. 


a PHYSICAL SYSTEM > \ i ` 
1 1 

! I 
I ! I I 
i 1 i 1 
I 1 i I 
I 1 i 1 
1 JOINT OF A ! 1 ! 
i WING SPAR ! i ! 
s » / 1 

U Ee ome TE n 3 il FEAE ET M 

I 1 I Y 
l i 7 t ! I ! 
I an "mid Magd Shell 1 i r 
LI \ 1 i 1 
I 1 I 
1 1 l I 
| | | 
LI 

r MODEL IDEALIZATION 4 E. MODEL IDEALIZATION " 

VSS T ae = uuu oe urere um iu - uu € um A MED Can Ew uw m ce aw 


Figure 1.5: Model Idealization of two case scenarios. 
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Figure 1.6: Examples relating the physical and mathematical models. 


For example, the aircraft structure, which has many components, must be separated 
into substructures in order to make the problem tractable: wings, fuselage, stabilizers, 
engines, landing gears, and so on. Each substructure is broken into its components: 
rings, ribs, spars, cover plates, actuators, etc. Eventually, those components become 
sufficiently simple in geometry and connectivity to allow a reasonable description through 
mathematical models based on the physics of the problem. It is possible to formulate a 
governing equation corresponding to a physical phenomenon and then obtain a system of 
equations through its discretization. This allows solving for the variables of interest such 
as deformations, stresses, and strains, as illustrated in Fig. 1.6. 


The finite element method can be applied to partial differential equations governing 
the physics of a problem or a purely mathematical problem with no relationship to 
a physical problem. In either case, a partial differential equation can be discretized, 
leading to a system of linear or nonlinear algebraic equations that, when solved, provide 
approximate solutions to the partial differential equation over the complicated domain. 


Figure 1.7 shows a heat transfer problem where the analytical solution is too cum- 


bersome to use even for a simple square domain, and how by applying FEA, the problem 
turns into a linear system of equations that can be easily solved by matrix inversion. 
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Figure 1.7: Benefits of FEM when compared to the analytical solution. 


1.2.4 Verification and Validation 


Figure 1.8 shows oversimplified steps of the FEM. The finite element method is an — (Felippa, 2007) 
approximate numerical method to solving partial differential equations for idealized 

physical systems, leading to unwanted errors. Many sources of errors could exist, and 

so these models must be adjusted using experimental measurements. The error in the 

solution is the amount by which the discrete solution fails to satisfy the discrete equations. 

There are two methods for error identification: verification and validation. The many (Anon. 20213) 
possible assumptions made throughout the finite element model could introduce errors. 
Unfortunately, these errors could be challenging to identify and can adversely affect the 

design process. 


IDEALIZATION DISCRETIZATION SOLUTION 


Physical 
System 


Modeling * Discretization * Solution error 


Figure 1.8: A simplified view of the physical simulation process to illustrate the finite 
element process. (Felippa, 2007) 
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Verification and validation, or known in its short form V&V, consists of developing 
computational models to make engineering predictions with quantified confidence. Verifi- 
cation and validation are a process of accumulating evidence of a finite element model’s 
correctness and accuracy. 


1. Validation refers to the process of checking whether the finite element simulations 
reflect real-world results. In other words, validation is the process of determining 
whether the finite element model precision represents the real system from the 
perspective of the projected uses of the model. 


2. Verification refers to the process of checking that the FEA was conducted properly. 
More precisely, verification is the process of defining that a model implementation 
correctly represents the developer's conceptual description of the model and its 
solution. 


It is common practice to verify finite element models and validate their analysis 
against experimental tests. Therefore, FEA must provide meaningful and valuable data 
against which to run experimental validations. Engineers need to verify the precision of 
the FEA models and validate the outcomes against experimental tests to guarantee a final 
model with the smallest possible error. The V&V processes reduce the time, cost, and risk 
associated with full-scale testing of the structures or system. 


Physical 
System 


TESTING 
SIMULATION ERROR 


Figure 1.9: Model updating process in the Physical FEM. (Felippa, 2007) 


Errors can occur in any of the steps in the finite element method: idealization, 
discretization, and solution. Consequently, it is essential to perform model verification 
and validation to identify errors and adjust to minimize these errors within the model. 
Model updating is a way of adjusting the finite element model to represent the physics 
better. The comparison of the discrete solution against the experiments helps find the 
parameters to adjust the FEM, as shown in Fig. 1.9. Since the minimization of errors 
is generally a nonlinear problem (even if the model is linear), the updating process is 
inherently an iterative one. Figure 1.10 shows an example of an aerospace composite 
structural component’s model validation process. 
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Figure 1.10: Example of model validation process of an aerospace composite structural 
component. 


1.2.5 Model Discretization 


The accuracy of the FEA depends largely on the choice of elements that best capture the 
physical system, the mesh (subdomains of the domain) type, and mesh density (number 
of subdomains per unit area). The development and improvement of many types of finite 
elements help to idealize the physical system's features. Each element formulation aims to 
approximate the physical response mathematically. Figure 1.11 shows several types of 
finite elements. Elements can be one-dimensional, two-dimensional, or three-dimensional. 
There are also unique elements with zero dimensionality, such as lumped springs or 
point masses. By using kinematic transformations, it is possible to expand the inherent 
dimensionality of the element. For example, it is possible to use a 1D element such as a bar, 
spar, or beam to build a model in 2D or 3D space. The solution within an element is either 
approximated with a linear polynomial (first-order element) or a quadratic polynomial 
(second-order element). 


These elements fall into several categories, and the focus of the discussion limits to 
structural analysis: 


1. 1D Bar Elements: These bar elements can have two nodes for linear approxi- 
mations and three nodes for quadratic approximations, and they aim to idealize 
three-dimensional structures that are long and slender and can only carry loads 
along the length of the bar. An example is a truss element. Bar elements can be 
oriented anywhere in 3D space. These elements transmit force along the length of 
the bar, and the nodes can only translate. 
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Figure 1.11: Typical finite element geometries in one through three dimensions. (Felippa, 
2007) 


2. 1D Beam Elements: Generally, there are two types of beam elements. These 
elements aim to idealize three-dimensional structures that mainly carry bending 
loads. One beam element type is intended for thick beams (Timoshenko beams) 
and can usually have either two nodes for linear approximation or three nodes for 
quadratic approximation. Slender beams are approximated with two nodes using a 
cubic interpolation function (Eutler-Bernoulli beam). Beams can also be oriented 
anywhere in 3D space. Each node of a beam element has six DOF elements allowing 
both translation and rotation. 


3. 2D Planar Elements: Planar elements can be linear triangular elements or second- 
order triangular elements and linear quadrilateral elements or second-order quadri- 
lateral elements. These planar elements aim to idealize three-dimensional structures 
that can be approximated using plane stress, plane strain, or axisymmetric assump- 
tions. These ideas will be further explored in this book. These nodes of these 
elements can only support translational DOF. 


4. 2D Planar Membrane Elements: Membrane elements can have any orientation in the 
3D space and help model three-dimensional structures that act as thin membranes 
(e.g., fabric). The elements can be triangular or quadrilaterals, just in the case of 2D 
planar elements. 


5. 2D Plate Elements: Plate elements are either triangular or quadrilateral elements, and 
it is possible to give them orientation in the 3D space. These elements can employ 
linear or quadratic approximations for the variable of interest. These elements are a 
suitable idealization of shell structures such as the walls of a pressure vessel and 
aircraft fuselages. These elements support translational and rotational DOFs. 


6. 3D Solid Brick Element: Brick, hex, or tetrahedral elements can be first-order or 
second-order, and the nodes only support translational DOF. These elements are 
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Figure 1.12: Example of typical idealization of primitive structural elements. (Felippa, 
2007) 


the most realistic representation of the solid object but can also be computationally 
expensive: 

(a) 3D Solid First-Order Tetrahedral (4 nodes) Element 

(b) 3D Solid First-Order Hex (6 nodes) Element 

(c) 3D Solid First-Order Brick (8 nodes) Element 

(d) 3D Solid Second-Order Tetrahedral (12 nodes) Element 

(e) 3D Solid Second-Order Brick (20 nodes) Element 


This book provides the elemental derivations for some of the above standard element 
formulations. The choice for the appropriate finite element depends on the idealization of 
the physical structural component. The degree of freedom at any point in the structural 
component influences the choice, as well. Figure 1.12 shows a few examples of mainly 
one-dimensional structural elements, and Fig. 1.13 shows a few examples of two- and 


three-dimensional structural elements. (Felippa, 2007) 
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Figure 1.13: Example of typical idealization of continuum elements. (Felippa, 2007) 


FEM for Solid Mechanics 


(1) Formulate systematic solution procedure for the PDE over a simple 
subdomain “elements” (element formulation) 


(2) Discretize the geometry into these subdomains (“pre-processing) 
E J T of elements 


D < di Trusses, beams 


ta 2 A » Plates, shells 
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Break-up CAD model into simple 
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based on physics (Mesh) | 


Figure 1.14: Mesh dependency on type of finite elements. 


Using a mesh with large element sizes, also known as using a coarse mesh, saves 
computational time, and it is a good initial step in ensuring that the model is working 
properly. In subsequent steps, mesh density is increased in critical areas. It is customary 
to compare the output results across various mesh refinements. One approach to make this 
comparison fruitful, is to plot the output for the variable of interest in the y-axis, and either 
the element size or the number of elements in the x-axis. After comparing a minimum of 
three successive solutions, the asymptotic behavior of the solution starts to emerge, and the 
changes in the solution between the various mesh densities become smaller. Figure 1.15 
shows an example of mesh refinement. This figure shows that as the number of elements 
increases, the stress concentration factor for a centered hole plate under axial tensile load 
approaches asymptotically to its theoretical value. 
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Figure 1.15: Approximating static stress concentration factor for an open whole plate 
under axial tensile load using Abaqus. 


There are two approaches employed with the purpose of increasing solution accuracy: 
h-refinement reduces the element size, and p-refinement increases the element order. While 
both approaches can increase computational time, computational efficiencies continue to 
increase as technology advances. Figure 1.16 shows an example of the mesh refinement 
for a global metric. In contrast, Fig. 1.17 shows an example of a mesh refinement. There 
is no need to refine the entire joint except in the bonded region, where damage is studied. 
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Figure 1.16: Idealizing the physical system to assist with FEA. 
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Figure 1.17: Example mesh refinement of an idealized composite joint and its solution. 
Courtesy of Jim Lua, Global Engineering, NJ. (Pham et al., 2021) 
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1.2.7 Basic FEM Steps 


There are three basic parts of the FEM, and each step is as important as the other: 


1. Preprocessor. In this step, the user prepares the geometry, cleaning it up, and mak- 
ing decisions on what details to include or omit. The geometry is then discretized 
using the most appropriate element types. Section properties are then assigned to 
the element types. For example, a beam may require the specification of a hollow 
circular cross-section. Material properties are specified and assigned to the regions 
of the models. Material models can include linear elastic isotropic or composites, 
plasticity, creep, etc. Loads and boundary conditions are specified according to the 
physics of the problem. Finally, the most appropriate solution type is selected: heat 
transfer, static analysis, dynamic analysis, etc. 


2. Processor. Includes all the steps in the finite element method (except for postpro- 
cessing). Element matrices based on the element formulation are determined. Next, 
the system of equations for the global system is assembled using local-to-global 
connectivity data by relating element nodes to the global nodes. Boundary condi- 
tions are imposed by identifying the restrained and unrestrained degrees of freedom. 
Lastly, using numerical techniques helps to solve algebraic equations for the primary 
variables at the nodes. 


3. Postprocessor. The solution to the system of equations is then utilized to post- 
process other variables of interest. Since each element uses linear or quadratic 
approximations, these polynomials can interpolate the solution within the element 
by using the nodal quantities found during the processing step. As a result, it is 
possible to find secondary variables, such as stress in solid mechanics problems or 
heat flux in a heat transfer problem. Lastly, the output data are processed for graphic 
display or printout. In the postprocessing step, the user performs verification and 
validation activities. Highly encouraged are hand-calculations, whenever possible, 
as a means to verify the numerical results. 


1.3 Boundary Value Problems 


The description of many physical problems such as heat transfer and solid mechanics 
consists of a boundary value problem described by partial differential equations and 
associated boundary conditions. These equations can be time-dependent or dependent, and 
can be linear or nonlinear in nature. A solution to a boundary value problem is a solution 
to a set of partial differential equations which must also satisfy boundary conditions. For 
example, the heat conduction equation in heat transfer must be solved for temperature 
while also considering that temperatures at the surface of the structure may be specified, 
so the solution must satisfy the problem boundary conditions and the domain’s governing 
equation. 
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1.3.1 Time Independent and Dependent Problems 


Physical problems are always time-dependent, but if the nature of the forcing function is 
slow relative to the response of the physical structure, then approximations can be made 
by considering the problem as time-independent. For example, a structure loaded very 
"slowly" can be considered time-independent and solved statically rather than considering 
time derivatives because, in this instance, the slow application of load is unlikely to 
cause a dynamic amplification. Specific to solid mechanics, type of time-independent 
analysis includes static and quasi-static. For the former, there is no need to include 
time explicitly. However, for quasi-static structural problems such as metal creep, rate- 
dependent plasticity or fatigue cycling, a realistic measure of time is required, but still, 
inertial forces are negligible. A steady-state analysis would be a time-independent analysis 
because it does not involve any time derivatives. 


Transient problems are problems where the primary variables of interest and the 
forcing functions depend on time and inertial forces. In these cases, the finite element 
formulation requires the calculation of time derivatives, and the equations of motion need 
to be integrated to march the solution over time. An example is a bar subjected to high 
heating, and the temperature response needs to be determined over time. Figure 1.18 
shows the analysis of bird-strike events to design turbo engine fan blades. (Goyal et al., 2013) 


BIRD-STRIKE DAMAGE ON ENGINE ARBITRARY EULER-LAGRANGIAN 


Figure 1.18: Modeling of Bird-Strike Events. (Goyal et al., 2013) 
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1.3.2 Linear and Nonlinear Problems 
(Anon., 2021b) 


A linear static analysis is one where a linear relationship exists between the applied forces 
and displacements. In structural problems, stresses must remain in the linear elastic range 
of the material, and the deformations are infinitesimally small so that linear geometry 
is sufficient. In linear static analysis, the stiffness matrix is constant, and therefore it is 
computationally efficient as the stiffness matrix needs to be inverted only once. 


A nonlinear analysis is one where the relationship between forces and displacements 
is nonlinear. The stiffness matrices in the nonlinear analysis can depend on the deflections, 
thus requiring a solution strategy capable of solving nonlinear problems, such as the 
Newton-Raphson procedure. Nonlinearity can result from: 


1. Geometric nonlinearity: Occurs in problems where large deformations occur. Many 
engineering applications such as metal forming, tire analysis, and medical device 
analysis require large deformation analysis. A few examples where small deforma- 
tion also requires geometric nonlinearity are cables, arches, and shells. 


2. Material nonlinearity: Involves the nonlinear behavior of a material based on a 
current deformation, deformation history, or rate of deformation. Examples of 
nonlinear material models are plasticity, elastoplasticity, and hyperelasticity (rubber 
and plastic materials). 


3. Follower loads: These are situations where the load depends on the deflection of 
a structure. For example, a load that remains tangent to the beam’s axis during 
deformation is a follower force because the load depends on the deflection and slope 
of the beam at that location where the load is applied. 


4. Contact: Usually requires nonlinear procedures, as deflections can depend upon 
other deflections within the same system. 


1.3.3 Applications to Structural Analysis 


For structural problems, it takes fifteen equations and fifteen unknowns to express the entire 

elasticity field. The reader is encouraged to review a complete treatment of elasticity from 

the previous volume Solution to Engineering Problems Using Computational Mechanics. (Goyal and Goyal, 2021) 
These fifteen unknowns are: 


1. Displacement field (3 unknowns):A displacement field d(x, y, z) is kinematically 
admissible if it is continuous and differentiable at all points in the domain. It is 
convenient to express the three unique displacements in vectorial form as follows 
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2. State of Strain (6 unknowns): A strain field E(x, y, z) is compatible if it is derived 
from a kinematically admissible displacement field through the strain-displacement 


relationships. 
Ex Ey Ey, 
E=| Ey Ey Ey (1.1) 
Ex, Ey, Ex 


It is convenient to express the six unique strain quantities in vectorial form as follows 


Exx €xx 
Eyy €yy 
EE Ew EN €zz 
E = = 
Eyz 2 eyz 
Ex; 2 €xz 
Exy 2 €xy 


3. State of Stress (6 unknowns): A stress field S(x, y, z) is statically admissible if it 
satisfies the equilibrium equations at all points in the domain. The stresses corre- 
sponding to small deformations and strains, in their three-dimensional Cartesian 
coordinate tensor notation, is 


Sx S xy Sxz Oxx Uy  Txz 
S= | Sy Sy Sy | | hy Oy Ty 
Sz S yz Szz Txz  Tyz Oz 


It is convenient to express the six stress quantities in vectorial form as follows 
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The fifteen equations governing the mechanics of a deformable body are as follows: 


1. Equilibrium equations (3 equations) must be satisfied at all point inside the body. 
The equations of equilibrium are the most fundamental equations and can be derived 
from integral form of Newton's law stating that the sum of all the forces acting on a 
differential element of the structure should vanish: 


d Oxx i OT | Ot 


dx ^ oy dz ps 
Ot Joy JTy are 
dx P dy az DE 
Ot, Ot, 00, n 
dx = dy dz apa 


where b,, b, and b; are the body forces per unit volume and x, y and Z the accelera- 


tions acting on the body. The components of this body force vector resolved in the ^ body force is a force that acts on the volume of a 
Cartesian coordinate system is body and can be defined as an external force acting 
throughout the body's mass. It has units of force 
by per unit volume and is constant within the element. 
b= b, Gravity is a body force, and the fictitious forces 
E b. associated with a non-inertial reference frame may 
z be viewed as body forces. Examples of common 
body forces include: (i) gravity, (ii) centrifugal 
force, (iii) coriolis force, (iv) acceleration felt when 
2. Strain-Displacement relationship (6 unknowns). The strain-displacement equations a vehicle changes velocity (as in braking), and (v) 
govern the deformation of the body at a point. The strain-displacement relationships — ? pressure gradient acting within a fluid (but not 
can be derived using differential geometry. For large deformations, the Green- Pressure alone). 
Lagrange strains are the appropriate strain measures, and the derivatives are defined 
relative to the undeformed configuration. The normal strain components in terms of 
the displacement gradients are follows See Appendix C. 


_au 1f/9UY (av a aw V 
x= Ox 2) Nox) "Vox EN 

| Ov 1f/0UY (av n aw \* 
ow By 2 oy oy oy 

.9w 1 f (au E av c aw V? 
“Oz AE dz dz 


For moderately small rotations, the angular distortion, or the shear strain about the 
x-axis can be approximated as 

_0V oW dU dU JVV JW OW 
Wor Oy 0r Oy 02 ay de 
dU JV dUdU dV AV AW OW 
hy = “a FEET ox" dy Ox 
_ dU JW JUAU 2O0VOV JW OW 
te = tr xz Ox dz) Ox à 


The units of the force vector are force per unit volume. 


https://gioumeh.com/product/engineering-problems-using-finite-element-methods-solution/ 


1.3 Boundary Value Problems 27 


For small deformation gradients, the equations reduce to their linear form: 


ðU oV JV 
ne pta 
_ov _dU JV 
ey = Dy f= cde m 
oW QU JW 


"wc jr Bac Toe 


3. Linear Stress-Strain relationship (6 unknowns) are known as constitutive relation- 
ships and relate stress to strain. These relationships can be found by performing 
laboratory testing. For elastic anisotropic materials the most general form is as 


follows: 
Sxx Du Di Dis Dis Dis Die exx 
Syy Di; Da Da Da Dos Dog yy 
S, L | Diz Doa D33 D34 D35 D36 Ez; 
Sy [| | Dis Doa D34 Dag Das Dag yz 
Sxz Dis Dos Dss Das Dss Dsg €xz 
Sxy Dig Dog Dag Das Ds6 Dee €xy 


Using virtual tests, it is possible to determine the relationship for an orthotropic 
material system in terms of strain-stress as follows: 


; WE  —Vay/Exe —ValEx 0 0 0 s 
i —Vyx/Eyy 1/Eyy  —Vyz/Eyy 0 0 0 S 
be Va /Eg  —Vey/Ez 1/Ez 0 0 0 S. 
>? 0 0 0 1/Gy 0 0 E 
2e B 0 0 0 0 1/G,, 0 Se 
x 0 0 0 0 0 Gy xy 


These engineering constants are determined experimentally. The elastic moduli Exx, 
Eyy, Ezz; shear moduli Gyy, G,;, and Gyz; and Poisson ratio nu,y, nuy;, nuy;. The 
compliance matrix above is a symmetric one. 


4. In the theory of elasticity, the boundary conditions on the surface of a volume are 
either specified displacements or specified tractions loads at a particular point on 
the surface. Cauchy's relationship needs to be used to relate the tractions on the 
surface to stresses inside the body using the following formula: t; = $ijn;, where 
i = x,y,z and n; are the components of a normal vector at a point on the surface of 
the specified traction. 


Hence, there are fifteen unknowns and fifteen highly coupled equations over a 
general domain. Although analytical solutions may be possible for simple geometries, in 
general, this is not the case. The finite element method consists of solving the elasticity 
field approximately by first solving for deflections, then using the deflections to solve 
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PRE-PROCESSING 


Static: Displacements, strain, stress 
Dynamic: Frequencies, Mode shapes 


Figure 1.19: Finite element analysis for a centered hole plate experiencing a pulling load. 


for strains, and then using the strain calculations to solve for stresses, thus solving all 
fifteen unknowns for the elastic field. However, the solution to the equivalent integral 
representation of the boundary value problem is the solution to the fifteen unknowns. 
Furthermore, the finite element method approximates the boundary value problem by 
solving the integral representation. Figure 1.19 shows how all these equations interact 
together in the FEM. 


1.3.4 Applications to Heat Transfer 


Partial differential equations in heat transfer problems are in terms of temperature. Heat 
source or sink, Q, over the volume is balanced by heat flow per unit area per unit time 
across surface q». 


I Qav - Í andA (1.2) 


The heat flow across the boundary surface is 


qx 
dn» —Q:n or q= qy (1.3) 
qz 


Using Fourier's Law, which states that the rate of flow of heat energy per unit area through 
a surface is proportional to the negative temperature gradient across the surface, it is 
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possible to obtain the following relationship: 


d 
dx 
Kxx kyy Kxz o Ax 
q= -DVT =- | ky ky ky 3. T= dy (1.4) 
Kxz YZ kzz y qz 
d 
dz 
Hence, 
d, —9:n—q'n-— -(DVT)'n (1.5) 


Now, use the Gauss theorem to convert the right-hand-side of the Eq. (??) into a volume 
integral: 


J QdV = ni qndA = ni Vl qav (1.6) 
V S=0V V 


Thus, the governing equation becomes: 
J| (c-vigav= [|| o+vrovrjav 00+ v'vr =0 (1.7) 
V V 


For orthotropic conductivity and extending the governing equation to a time-dependent 
problem as follows: 


oT 
Ox 
a d ə kx 0 0 a 
vipvr+0={ 5 5 x] 0 ky 0 — p+0= 
m ee 2s 0 0 kz 2 (1.8) 
T 
k aT k OT oa E ql : 
XX EN, yy ay ZZ a = Pcp ðt 


where p represents the density, cp the specific heat capacity of the material, T temperature, 
and t time. For isotropic materials only, the conductivity k is specified as follows: 


oT T APT oT 
EE ay? | EE += Peps, (1.9) 


where the associated boundary conditions are: 


1. Prescribed surface temperature: The surface temperature is prescribed as a function 
of position and time 


T(x,y,z,t) = T (x,y,z,t), 


where T(x,y,z,t) is a prescribed temperature function of the position and the time. 


https://gioumeh.com/product/engineering-problems-using-finite-element-methods-solution/ 


30 Chapter 1. Introduction to Finite Element Method 


2. Adiabatic surface: When there is no heat flux through the boundary surface, the 
boundary conditions are as follows: 


ITY) _ y 
on 7 
3. Heat flux specified at the surface: When the heat flux is prescribed on the boundary 
surface, the boundary conditions are: 


OT (x, y,z,t _ 
—k E ) = qa(x,y,z,t) 
n 


where Gn (x, y, z, t) denotes the heat flux prescribed on the boundary surface. 


4. Heat transfer by convection: In many cases, the surface of the solid is in contact 
with fluid or gas. When the heat transfer between the boundary surface of the solid 
and the surrounding medium occurs by convection, it results in the following mixed 
boundary condition: 


oT 
—k—— +5 = he (T — Too 
a5 Qs c ( ) 
where qs represents heat generation on the surface, T the surface temperature, To. 
the temperature of the fluid or gas surrounding the structure, and h, the heat transfer 
coefficient. 


5. Heat transfer by radiation: A body at absolute temperature T can exchange heat by 
radiation. 


oT 
—k-— t4; = 0€ (T^ — T^ 
an ds ( ) 
where o represents the Stefan-Boltzmann constant, e the emissivity of the surface, 
and 7; the temperature of the medium. Since this equation is nonlinear, it may be 
convenient to derive a linear approximation. If the temperature difference between 
T and T. is small, then the boundary condition simplifies to: 


oT 
HA +qs = 408 T2 (T — Ta) 


6. Heat transfer by convection and radiation: When both radiation and convection are 
considered, the boundary condition becomes: 


oT 
kor +4s = he (T — Te) + oe (T! — T-*) 


Time-dependent problems require specified initial conditions are required and tem- 
perature distribution in the solid at the initial time. If the temperature reaches steady-state, 
it is not necessary to specify an initial condition. 
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